#!/usr/bin/env python2
from pylab import * 

name = 'test'
#name = 'test2'
#name = 'MLE'


cache_covariance = 0
cache_covariance_det = 0
cache_covariance_inv = 0

def normal_distribution_pdf(point, mean, covariance):
    global cache_covariance
    global cache_covariance_det
    global cache_covariance_inv
    if  cache_covariance is not covariance:
        cache_covariance_det = det(covariance)
        cache_covariance_inv = inv(covariance)
        cache_covariance = covariance
    point_offset = point - mean
    tmp0 = 1 / (pow(2 * pi, point.size / 2.0) * sqrt(abs(cache_covariance_det)))
    tmp1 = dot(dot(point_offset, cache_covariance_inv), point_offset)
    return tmp0 * exp(-0.5 * tmp1)


def load_normal_distribution(file_name):
    mean = loadtxt(file_name + '_mean')
    covariance = loadtxt(file_name + '_covariance')
    return mean, covariance

def plot_normal_distribution(mean, covariance, grid_):
    def tmpfunc(point):
        return normal_distribution_pdf(point, mean, covariance)
    return apply_along_axis(tmpfunc, 2, grid_)


N = 50
r = 4.0
x = linspace(-r, r, N)
y = linspace(-r, r, N)

X, Y = meshgrid(x, y)
grid_ = dstack((X, Y))

#############
mean, covariance = load_normal_distribution(name + '_shooting_distribution0')
data0 = plot_normal_distribution(mean, covariance, grid_)
mean, covariance = load_normal_distribution(name + '_shooting_distribution1')
data1 = plot_normal_distribution(mean, covariance, grid_)
#pcolor(X, Y, data1 / (data0 + data1))
pcolor(X, Y, data0 + data1)
######################

shots1 = loadtxt('test_shots1')
plot(shots1[:, 0], shots1[:, 1], 'wx')
shots2 = loadtxt('test_shots2')
plot(shots2[:, 0], shots2[:, 1], 'r+')
colorbar()

show()  
